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This paper describes two methods of trajectory optimization to obtain an optimal trajectory of minimum- 
fuel-to-climb for an aircraft. The first method is based on the adjoint method, and the second method is 
based on a direct trajectory optimization method using a Chebyshev polynomial approximation and cubic 
spine approximation. The approximate optimal trajectory will be compared with the adjoint-based optimal 
trajectory which is considered as the true optimal solution of the trajectory optimization problem. The adjoint- 
based optimization problem leads to a singular optimal control solution which results in a bang-singular-bang 
optimal control. 


I. Introduction 

With the increase in fuel cost in the aviation industry, focus on energy efficiency becomes an important driver in 
aircraft design and operation. Lightweight aircraft design is a current trend in the aviation industry that addresses the 
need for increased energy effficiency in aviation systems. Improved propulsion systems also contribute significantly to 
the goal of energy efficiency. Another way of improving energy efficiency is by means of strategic operational design 
of mission profiles to achieve reduced energy consumption. Strategic mission performance addresses optimization 
of flight path that achieve a certain mission objective. In this paper, we are interested in exploring strategic mission 
design for minimizing energy consumption during climb. 

Trajectory optimization for climb has been well-studied. 3 The traditional adjoint-based trajectory optimization 
method can be used to obtain an optimal trajectory. This method requires solving a two-point boundary-value problem 
(2PBVP) that enforces initial conditions in the state equations and terminal conditions in the adjoint equations. The 
adjoint-based method is considered an indirect optimization method which provides a correct optimal solution. The 
solution of a 2PBVP can sometimes be difficult and the development of the adjoint equations can be tedious. Thus, 
various approximate adjoint-based optimization methods as well as direct optimization methods have been developed. 
Adaptive critic is an approximate dynamic programming method that can be used to solve a 2PBVP. For direct op- 
timization, the direct collocation method has gained popularity in recent years and has been used in many trajectory 
optimization problems. Approximate direct optimization methods using cubic splines have also been developed. In 
this paper, we conduct three methods of trajectory optimization to obtain an optimal trajectory of minimum-fuel- 
to-climb for an aircraft. The first method is based on the adjoint method. The second method is based on a direct 
trajectory optimization method using a Chebyshev polynomial approximation. The third method is based on a direct 
trajectory optimization using a cubic spline approximation. The approximate optimal trajectory will be compared 
with the adjoint-based optimal trajectory which is considered as the true optimal solution of the trajectory optimiza- 
tion problem. The adjoint-based optimization problem leads to a singular optimal control solution which results in a 
bang-singular-bang optimal control. 4 
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II. Adjoint-Based Trajectory Optimization for Minimum Fuel-to-Climb 


Consider the standard aircraft point-mass model described by the equations of motion in a vertical plane with the 
assumption of the thrust angle aligned with the aircraft x-axis 


V = 


7 = 


T — D — W sin y 
m 

L — W cos y 


mV 
h = V sin y 
W = —cT 

subject to the initial conditions V ( t , ) = V;, y(f,-) = y, /z (/,-) = /z,-, and W (t,-) = W) where f; is the initial time. 
Consider a climb trajectory with minimum fuel consumption. The cost functional to be minimized is 


J = 


rtf 


cTdt 


( 1 ) 

( 2 ) 

(3) 

(4) 


(5) 


where c is the thrust specific fuel consumption (TSFC) , T is the engine thrust, and tr is free terminal time at which 
the aircraft reaches a specified airspeed V (tf) = V f and altitude h ( tf ) = /z/. 

The TSFC and engine thrust are generally functions of airspeed and altitude. 

Consider a simplified, quasi-steady flight path model by assuming thaty ss 0 and sin y sa y for small y. Defining F 
as the specific excess thrust (SET) 


F(V,h,W) 


T (V,h) —D(V,h) 
W 


(6) 


the equations of motion can be re-expressed as 


v = g (F-r) 
ii = Vy 
W = -cT 


(7) 

( 8 ) 

(9) 


We formulate the Hamiltonian function as 


H = cT+ X v g (F — y)+XhVy— X w cT 

The adjoint equations are obtained as 

• dH „ ,d(cT) , dF , 

Av = = (A w - 1) -dy--*vgjy-kr 

j dH d(cT) dF 

■ dH dF F 

Xw = ~dW = ~ Xv8 dW =lv8 W 

subject to the transversality conditions at the terminal time Aw (tf) = 0 and 

cT + Xvg(F-y)+X h Vy\ t=tf = 0 
The optimal control is determined from 


dH 

dy 


= —Xyg + XhV 


( 10 ) 

( 11 ) 

( 12 ) 

(13) 

(14) 

(15) 


Note that we cannot solve for y from dH /dy because y is not present. This is a special case of singular optimal 
control which is defined when the control variable appears linearly in the Hamiltonian as seen in Eq. (10). Let us 
rewrite the Hamiltonian as follows: 

H = (l-X w )cT + X v gF + Sy (16) 
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where S is called a switching function which in this case is 


S = = — Xyg + A/,V (17) 

To minimize the Hamiltonian, we invoke the Pontryagin’s maximum principle 

y* = arg 7 ey ad min H (18) 

Thus, the optimal control is one that minimizes the Hamiltonian. This results in the following optimal control 

! Ymax 5 ’ 0 

Ysing S = 0 (19) 

Ymin S 0 

This type of control is called a bang-singular-bang control, which implies that the control switches between two 
extremal values depending on the switching function S. A natural question would arise as to what the optimal control 
y would be if S = 0. The optimal control that corresponds to S = 0 is called a singular optimal control, which can 
be determined by taking successive time derivatives of S and setting them to zero until the control variable appears. 
Taking the first time derivative of S yields 


S = — Xyg-\-Xi,V + A/,V = 0 

Substituting in the adjoint equations Xy and Xi, yields 

.. , , <9 (cT) . o dF d (cT ) , dF 

~g{Xw - 1) w - 1) — ^ = ® 


Since 5 = 0, therefore 

Substituting A/, into Eq. (21) yields 

~ 2 / dF V 2 dF 


Xh = 


Xyg 

V 


g 

r 

v 


F + V 

dV g dh 


+ (Xw ~ 1 ) 


d(cT) d(cT ) 
dh 1 dV 


= 0 


Since tf is free, the Hamiltonian is zero at all times. Therefore 

XygF 


Xw 1 


cT 


( 20 ) 


( 21 ) 


( 22 ) 


(23) 


(24) 


Upon substitution, the singular arc in the state space along which the singular optimal control is applied is obtained 


as 


, . dF V 2 dF FV 

mh) = F + v^-- Th -- 


d (cT) V d ( cT ) 


dV 


g dh 


= 0 


(25) 


It is of interest to note that the first three terms of f(V,h) turn out to be the singular arc for a minimum-time-to- 
climb trajectory optimization problem. Thus, the minimum-fuel-to-climb trajectory can be thought of as the trajectory 
that minimizes both the time to climb as well as the fuel during climb. 

The singular arc is a unique trajectory in the V — h plane. Differentiating the curve yields 


df d f 

2 dV + -Frdh = 0 
dV dh 

The singular optimal control, which is the flight path angle y, is obtained from 


dV 

dh 


V _ g(F-y) 
h Vy 


(26) 


( 27 ) 
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which yields 


Ysing 


V dV 

~~n~Y = F 
g an 



df 

dh 


(28) 


For the initial segment, an initial value of y = y max is guessed. From the initial conditions, V and h are integrated 
forward in time until the trajectory intersects the singular arc at a flight path angle y = Ysing- If Ymax / Ysing - then a new 
value of y = y max is chosen with y max = Ysing- The solution is iterated until y max = Ysing- 

Similarly, for the final segment, a final value of y = y„„„ is guessed on the singular arc. From the initial values of V 
and h on the singular arc, V and h are integrated forward in time until h ( tf ) = hf. If V (tf) ^ Vj, then a new value of 
y = Ymin is chosen on the singular arc such that V new = Void + V/ — V (tf) and the solution is iterated until V (tf) = V/ 
and h (tf) = hf. 

Since the aircraft’s aerodynamic coefficients as well as the engine thrust and thrust specific fuel consumption are 
dependent on the Mach number M = V /a(h) where a is the speed of sound, it is expedient to express the singular arc 
as a function of the Mach number and altitude as 


f(M,h)=F + 


dF 

M- h 

dM 


M 2 a da 

8 dh 


dF 

~dM 


M 2 a 2 dF 

8 dh 


F / M 2 ada\ d(cT) 
cT \ + g dh) dM 


M 4 2 a 2 d (cT) 

8 dh 


(29) 


III. Approximate Trajectory Optimization Using Chebyshev Polynomials 

The optimal trajectory is a unique function in the V ~h plane. Given that the initial and terminal values of V and 
h are specified, a candidate function can be constructed to approximate an optimal trajectory. Toward that end, we 
propose the use of Chebyshev for approximating a trajectory optimization problem. 


A. Chebyshev Polynomial Approximation 

Any sufficiently smooth function f (xj £ r t? can be expanded as a Taylor’s series about some x = xo 

f(x)=f {xo) + V/* (x 0 ) (x - x 0 ) + * (x - x 0 ) T V 2 /* (x 0 ) (x - x 0 ) + ■ ■ ■ (30) 

Then / (x) can be represented as 


f(x) = ©* T 4>(x) + e(x) (31) 

where ©* £ K"' x W is a matrix of constant but unknown coefficients, <f> (x) £ R 9 is a vector of regressors in terms of 
monomials of x 


<F(x) = 


X\ x 2 


Xp x\ X\X 2 


q— 1 
X\X^ 



(32) 


and e (x) is a function approximation error or truncation error which depends on x. 
f (x) is then approximated by a polynomial of q ,h degree 


f{x)=Pq(x)=& ( t>{x) 


(33) 


where 0 £ R"' x R” is the estimate of 0*. 

The coefficients © can be computed using various least-squares methods such as the batch least-squares, least- 
squares gradient method, or recursive least-squares (RLS) method. Note that since 0 <J>(x) is an approximation of an 
unknown function / (x), the approximation error will not be asymptotic regardless whether or not <J> (x) is a persistently 
exciting (PE) signal. 

The Weierstrass theorem 1 states that given any sufficiently smooth function / (x) £ ^ [a, b] and So > 0, there exist 
a polynomial p q (x) for some sufficiently large q such that ||/(x) — < £o- This means that any sufficiently 

smooth function can be approximated by a polynomial of q ,h degree. Then the function approximation error could be 
made sufficiently small on a compact domain of x such that sup^g^ ||e (x) || < £o for all x £ $> C R". 

There are several types of polynomial approximation of a function. The regular polynomial regression using 
monomials as basis functions is frequently used for function approximation. Orthogonality is a property of a function 
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that belongs to a metric space endowed with an inner product. Given two functions g (x) and h (x), then g (x) and h (x) 
are orthogonal to each other if their inner product is zero. That is 


(g(x),h(x)) =0 


(34) 


where (.,.) is an inner product operator that takes on two real-valued functions and returns a constant. The inner 
product has a concrete mathematical definition depending on the class of functions and the type of metric space. 

A regular polynomial does not possess this orthogonality property. In contrast, certain classes of polynomials 
such as Chebyshev polynomials and Legendre polynomials are known as orthogonal polynomials. One advantage of 
an orthogonal polynomial over a regular polynomial is that an orthogonal polynomial can provide a better function 
approximation than a regular polynomial of the same degree. Because the Chebyshev orthogonal polynomials are 
considered as being nearly optimal polynomials 1 that can accurately approximate functions, neural networks based on 
the Chebyshev orthogonal polynomials have been developed. 2 

The Chebyshev orthogonal polynomials are special polynomial functions that are associated with the solution of a 
class of Sturm-Liouville differential equations 


(1-x 2 ) 


dry dy 
dx 2 dx 


+ n 2 y = 0 


(35) 


forx € [—1, 1]. 

This differential equation is known as the Chebyshev differential equation of the first kind. The Chebyshev differ- 
ential equation of the second kind is of the form 


( 1 ~ x 2 ) ^- 3 x t +n{n+2)y = ° 


(36) 


The Chebyshev orthogonal polynomials of the first kind are given by a generating function 


Tn + 1 (x) = 2xT n (x) - 1 (x) 


(37) 


where Tq (x) = 1 and 7) (x) = x. 

The first several terms of the Chebyshev orthogonal polynomials are given as follows: 


To(x) = 1 
T\ (x) = x 
Ti (x) = 2x 2 — 1 
7) (x) = 4x 3 — 3x 
7) (x) = 8x 4 — 8x 2 + 1 


(38) 


A solution of a Sturm-Liouville differential equation constitutes an orthogonal basis that spans any Hilbert space 
which is a complete, inner product space with a weighted inner product definition 


(g (x) ,h(x)) = w{x)g (x) h (x) dx 


(39) 


The Chebyshev polynomials are orthogonal with respect to a weighting function 

1 


'(*) = 


Vl — x 2 


such that 


(T„ (x) , T m (x)) = 


" l T n (x) T m (x) 
-1 Vl — x 2 


0 n ^ m 
dx= l K n=m = 0 


(40) 


(41) 


Any subspace -9 in an inner product space 9 has an orthogonal complement S? 1 such that their bases completely 
span the inner product space 9 . Therefore 

9 = (42) 
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Since the Chebyshev polynomials are orthogonal, they form a complete basis for an real-valued function. This 
implies that any unknown, non-singular function f (x) can be approximated by a Chebyshev polynomial of degree n. 

n 

f (*) = 0o7o (*) + 0i T\ (x) H h e n T n (x) = Y, OiTi ( x ) = © T< £ (*) (43) 

i= 1 


where the polynomial coefficients are approximated by 


0, 



f ( x )Tj (x) 

Vi — x 2 


dx 





(44) 


forx € [—1, 1], 

For higher dimensions when x G R", multi-dimensional Chebyshev orthogonal polynomials can be obtained. 2 


B. Trajectory Optimization with Chebyshev Polynomials 

Suppose a candidate trajectory in the V — h plane can be approximated by Chebyshev polynomials as 

V = ® T <t>(h) 


(45) 


where © is a vector of the unknown coefficients of the Chebyshev polynomials <t> (/;). 

The Chebyshev polynomials of unknown coefficients are subject to the boundary conditions at the initial and 
terminal times V (hi) = V l and V (hf) = V f. 

The trajectory optimization problem now becomes a direct search problem of finding © that minimizes a cost 
function. Toward that end, the cost function in t must be modified into a cost function in h as follows: 


rtf 


ftf 


n f cT 


J= cTdt = — / Wdt= / dh 

J tj J t : J hj h 


(46) 


where h is expressed as 


Therefore, the cost function is a function of © 


h = 


FV 


1 4. VdV 
1 + Rdh 


"f cT 


n®)= FT, 1 + --7T dh 


I hj 


FV 


V dV 


8 dh 


subject to equality constraints 


v l = e T <v(h l ) 

Vf = © T< f> ( hf ) 


(47) 


(48) 


(49) 

(50) 


Note that the integrand is a function of the unknown parameter ©. Therefore, the optimization is tantamount to 
computing the gradient of J with respect to ©. Thus 


dJ d T h f cT ( VdV\ JI 

dW~dWj hi Fv{ + -glh) 

= _ [ hf W_l F + V d_l +y d_Fd_h _FV 
J hi (FV) 2 1 dV dhdV cT 


d (cT) d (cT) dh 
dV + dh dV 


\ dV 

J <9© T 


, V dV\ „ 

1 H — ) dh 

8 dh 


h f cT d 


I hi Fvde T \ g dh 


v - d Zuh 


(51) 


A steepest descent method for convex optimization can be implemented to search for the optimal © that minimizes 
J as follows: 

©,-+!=©,• -6^(0,-) (52) 
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IV. Numerical Implementation 


The minimum-fuel-to-climb trajectory optimization problem is implemented using a Generic Transport Model 
(GTM) with the following thrust and TSFC characteristics: 


Tmax(M,h) = (— 77,03 1M+ 100,986) 5 (A) ^1 + ^M 2 ^) 7 ~‘ (53) 

c{MJi) = (0.45642M + 0.26156) ^9 (h) (1 +0.2M 2 ) (54) 

where M is the Mach number, 5 = p/psl is the pressure ratio, and 0 = T /T$l is the temperature ratio at altitude to 
sea level, with psl = 21 16.2 psf and Tsl = 518.69 °R are the standard sea level values. 

The aircraft initial weight is 200,000 lbs with the initial climb-out airspeed at Mach 0.2 at sea level and the terminal 


airspeed at Mach 0.8 at 35,000 ft. The clean-wing aerodynamics are described by 

C L {cc,M)=C Lo (M)+C La {M)a (55) 

C D (a,M) = C Dq (M)+K(M)C 2 {a,M) (56) 

where 

C Lo (M) = 0.27374M 3 - 0.24052M 2 + 0.095875 M + 0.034387 (57) 

C La (M) = 4.6340 M 3 - 3.1363 M 2 + 1.5094M + 4.5756 (58) 

C Dq (M) = 0.016454M 3 - 0.012294M 2 + 0.005823 \M + 0.015586 (59) 

K(M) = -0.019918M 3 -0.0035320M 2 -0.0064967M + 0.065037 (60) 

Below 1000 ft, high-lift flap extensions are deployed. 


The singular arc is obtained by solving Eq. (29). The other two climb segments from Mach 0.2 at sea level to the 
singular arc and from the singular arc to Mach 0.8 at 35,000 ft are obtained by trial-and-error. The two climb segments 
are at a constant flight path angle of 9.8035° for the first segment from sea level and 9.8035° for the third segment 
to 35,000 ft. The optimal trajectory and flight path angle are shown in Figs. 1 and 2. This climb results in a fuel 
consumption of 3,067 lb for a time to climb of 11.02 min. Figures 4 to 7 show the time histories of the altitude, Mach 
number, flight path angle, and aircraft weight change. 




Fig. 1 - Optimal Minimum-Fuel Trajectory Fig. 2 - Optimal Minimum-Fuel Flight Path Angle 
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Fig. 3 - Optimal Minimum-Fuel Aircraft Weight Change 




Fig. 4 - Altitude 


Fig. 5 - Mach Number 




Fig. 6 - Flight Path Angle Fig. 7 - Aircraft Weight 

To conduct the approximate trajectory optimization using the Chebyshev polynomials, the altitude is normalized 
by a variable x £ [—1,1]. The initial and terminal constraints on the coefficients of the first term and last term of the 
Chebyshev polynomials are imposed, respectively. These two constraints allows the coefficients of the first term and 
last term to be solved in terms of the remaining coefficients. The solution is feasible only for even-degree Chebyshev 
polynomials since odd-degee Chebyshev polynomials results in an inconsistency in the solution. 
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The approximate optimization is conducted using a random search method with the coefficients randomized to 
be between -0. 1 and 0. 1 . This range appears to provide better solutions than a larger range. Several degrees of the 
Chebyshev polynomials are used. Each degree of the Chebyshev polynomials is randomized by 10,000 iterations. At 
each iteration, any candidate trajectory that exceeds the Mach number range is discarded. The retained trajectories are 
used to compute the minimum fuel solution. 

Figures 8 to 10 are the plots of the approximate optimal trajectories using the Chebyshev polynomials of degrees 
4, 6, 8, and 10. The best solution is that corresponding to the 8-degree Chebyshev polynomial. The minimum fuel 
burns are given in Table 1 . 




Fig. 8 - Optimal Minimum-Fuel Trajectory 


Fig. 9 - Optimal Minimum-Fuel Flight Path Angle 



Fig. 10 - Optimal Minimum-Fuel Aircraft Weight Change 


Chebyshev Polynomial Degree 

Minimum Fuel Burn, lbs 

4 

3,640 

6 

3,284 

8 

3,143 

10 

3,232 


Table 1 - Minimum Fuel Burn with Approximate Optimal Trajectory by Chebyshev Polynomials 

While the fuel burn seems to be reasonably minimized by the Chebyshev polynomials, the approximate trajectories 
are highly oscillatory typical with polynomial approximation. From a practical consideration, this may not be entirely 
desirable. Thus, this approach does not seem to offer a significant benefit. 
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Another approach that could improve the trajectories is the use of piecewise continuous Chebyshev polynomials 
of low order to prevent oscillations. This would be effectively equivalent to using a cubic spline polynomial approxi- 
mation. This approach is implemented to see if the trajectory optimization could be improved. 


V. Trajectory Optimization by Cubic Spline Approximation Using Evolutionary Algorithm 

In this work the problem of finding the minimum-fuel-to-climb has been cast in such a way that an optimal solution 
can be found and in such cases this calculation is the approach to take. Yet for more complex setups - for example 
taking into account variable winds or aircraft flexibility - it is not feasible to calculate optimal solutions. For these 
situations using an optimization algorithm such as an Evolutionary Algorithm (EA) can be the best approach. 1 Here 
we describe our implementation of an EA-based approach for optimizing the minimum-fuel-to-climb. 

Optimizing the trajectory requires defining the space of trajectories to be optimized by creating a representation 
of a trajectory and constructing an objective function against which to optimize. Here the basic objective function is 
to minimize the fuel used during 30 minutes of flight time in going from a start state of an altitude of 100 ft and a 
speed of Mach 0.2 to a cruise altitude of 35,000 ft and a speed of Mach 0.8. To achieve this, the aircraft is set to use 
maximum throttle and the one control parameter available is that of setting the flight path angle throughout the flight. 
Flight path angle values are constrained to fall in the range [0°, 10°] and are sampled at a rate of 10 samples/s. Thus a 
candidate solution consists of a sequence of 18001 values in the range [0°, 10°]. 

From initial runs we noticed that what typically happens is that an aircraft reaches 35000 ft and a speed of less 
than Mach 0.8 and the unoptimized flight path angle had it continuing to gain altitude. Since creating a sequence of 
flight path angle values that correctly levels off at 35,000 ft until Mach 0.8 is reached adds significant difficulty to the 
problem, to make it easier we added a kind of autopilot. This autopilot is engaged when a candidate solution guides 
an aircraft to 35,000 ft, at which point the autopilot is engaged and all future flight path angle values are set to 0. 

The basic objective function used is to minimize the fuel used and these are measured in pounds. Since an EA 
starts with randomly generated candidate solutions, and these are not expected to achieve the target end state, the ob- 
jective function must include ways of scoring candidate solutions which fail to achieve the target end state. Candidate 
solutions which did not reach the desired altitude had their score penalized by the number of feet they differed from 
35000. Those which failed to reach the target speed were penalized by 10000 times how far their speed differed from 
the target cruise speed of Mach 0.8. Finally, to encourage resulting trajectories to be monotonically increasing in 
speed, a candidate solution’s score was given a SpeedPenalty of 1 for every time step in which it was slower than its 
previously-reached maximum speed. This objective function is summarized as follows: 

Score = Fuel Use - F |35000 — Final Altitude \ + 10000 x |0.8 — Final Mach\ + Speed Penalty (61) 

To represent sequences of flight path angle values we use a sequence of one or more cubic Bezier splines. Bezier 
splines are a type of smooth, parametric curve which is commonly used in creating trajectories, and by constraining 
the knots to fall in the range[0°, 10°], this results in flight path angle values being limited to the range [0°, 10°]. A 
cubic Bezier curve is specified by four knots (or control points) and a spline is just a sequence of these curves with 
an additional knot needed for each additional curve segment. While a single spline will produce a smooth curve, 
trajectories during climb often have up to three stages, with discontinuity between stages. To be able to create such 
trajectories, our representation can be configured to use a pre-specified number of splines connected in sequence, 
with a parameter for each spline indicating what percentage of the trajectory it specifies. The resulting representation 
system allows the user to specify at run-time both the number of splines to use and the number of knots per spline. 

This spline -based encoding of flight path angle values results in a vector of real values equal to the number of 
splines times the number of knots per spline. The EA we use for optimizing these values is the Age-Layered Population 
Structure (ALPS) EA, 1 in which multiple instances of a traditional EA are run in a sequence of layers. Individuals in 
the population have an age which is determined by the number of generations since their oldest ancestor was created 
from random. Each layer has a maximum age and at regular intervals the first layer is re-initialized with a new, 
randomly generated population. 

We use 12 layers of 40 individuals with an elitism of 5 in each layer. The maximum age for each layer was assigned 
based on a fibonacci-based aging scheme which have found works well and the maximum ages for each layer of were: 
3, 6, 9, 15, 24, 39, 63, 102, 165, 267 and 432. New individuals are created using mutation or recombination, with 
equal probability. In selecting parents, the first parent is chosen through a tournament selection with a tournament size 
of 5 and the second parent (only needed for recombination) is chosen at random. With mutation, either 1-4 genes are 
randomly selected to be mutated or all genes are mutated. Gene values are mutated by selecting a mutation size from 
1% to 0.000001% (by randomly selecting one of x 1 0 2 , x 1 0 3 , . . ., x 10 6 ) of the difference between the minimum 

10 of 12 


American Institute of Aeronautics and Astronautics 



and maximum values for that gene and then this is used as the variance for a random number with a normal distribution 
which is the mutation amount for that gene. For recombination, a value for each gene i is selected at random from the 
range [Pp, + (P\ , — Po,,) , Pi.i\, using a uniform distribution. 


Fuel Usage: 1 and 2 Splines 
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Fig. 1 1 - Fuel Usage for 1 and 2 Splines 


gamma vs time 



Fig. 13 - Flight Path Angle vs. Time 


Mach vs altitude 



Fig. 15 - Mach Number 


Fuel Usage: 3 and 4 Splines 



gamma vs altitude 



Fig. 14 - Flight Path Angle 


alpha vs altitude 



alpha 


Fig. 16 - Angle of Attack 


To find near-optimal trajectories we performed 5 runs using each combination of 1, 2, 3 and 4 splines with 6, 8, 10 
and 12 knots/spline. A subset of our results are presented in Figs. 1 1 and 12. These results show that poor results are 
achieved using only a single spline, or if only 6 knots/spline are used. All combinations of 2, 3 and 4 slines with 8 or 
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more knots performed fairly similarly and use just under 4259 lbs of fuel. One run with 4 splines and 12 knots came 
up with a slightly different solution which uses about 4255 lbs of fuel. All of these best solutions reach cruise altitude 
and speed in just under 20 minutes. 

The behavior of the best evolved trajectories are shown in the graphs in Figs. 13 and 14. The graph in Figure 13 
shows the flight path angle values for the best evolved trajectories using 1, 2, 3 and 4 splines. Most trajectories with 2, 
3 and 4 splines have almost identical flight path angle values up until around time step 8000. At this point, the aircraft 
have achieved the target cruise altitude and the remaining flight path angle values are ignore. Thus these trajectories 
are effectively identical, just as their scores on the objective function are nearly identical. Similarly, the behavior of 
these evolved trajectories is almost identical in looking at their flight path angle versus altitude, Mach versus altitude 
and angle of attack versus altitude curves (Figs. 15 and 16). 

The two exceptions to the results described above are the results for a single spline and for one run with 4 splines 
and 12 knots. All runs with a single spline were noticeably worse than the other runs, with adding more knots resulting 
in improved performance. Interestingly, the best solution found with 4 splines differs from the other solutions by 
starting with a small flight path angle value and rapidly increasing to near maximum, whereas the other solutions all 
start near the maximum. In addition, after it has reduced its flight path angle value it starts to bring it up about a minute 
after the other best solutions. This resulted in using 4255 lbs of fuel, which is about 4 lbs less than the other solutions. 
The fuel burn using cubic spline approximation is higher than the Chebyshev polynomial approximation due to the 
fact that Mach 0.8 is not reached when the final altitude is reached. This necessitates a constant altitude cruise until 
Mach 0.8 is reached. This constant altitude cruise results in a fuel burn penalty. Future work will examine techniques 
that could satisfy both terminal time constraints simultaneously. 


VI. Conclusion 

This paper presents a minimum-fuel-to-climb trajectory optimization using three methods. The first method is 
an “exact” method based on an optimal control formulation using the adjoint method. A singular control is obtained 
from this formulation. The second method is based on a Chebyshev polynomial approximation and the third method 
is based on a cubic spline approximation. The minimum fuel burn obtained by the optimal control method is “true” 
minimum, whereas the fuel burns from the two approximate methods are greater than the optimal value of the fuel 
burn. The Chebyshev polynomial method can provide a trajectory that could achieve a fuel burn near the optimal 
value, but it suffers from a problem with oscillatory trajectories. The spline method provides smoother trajectories 
than the Chebyshev polynomial method, but it cannot mutually meet the terminal time constraints, thus resulting in a 
higher fuel consumption. 
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